function [J] = jacobian(func, x, vars)
%JACOBIAN - Calculate the Jacobian matrix of a function.
%
% Inputs:
%   func - the function whose Jacobian we want to calculate.
%   x - the point at which to calculate the Jacobian.
%   vars - 
%
% Outputs:
%   J - the Jacobian matrix of func at x.
%
    % calculate the Jacobian
    % the column number of
    if nargin < 2
        error("Not enough input arguments!");
    end
    if nargin == 3
        if max(vars) > length(x)
            error("The function does not have so many arguments!");
        end
        J = [];
        for nvar = vars
            h = zeros(size(x));
            h(nvar) = sqrt(eps(1));
            J = [J, (func(x + h) - func(x - h)) / (2 * h(nvar))];
            % J = [J, (-func(x+2*h) + 8*func(x+h) - 8*func(x-h) + func(x-2*h)) / (12 * h(nvar))];
        end
    else
        J = jacobian(func, x, 1:length(x));
    end
end
